﻿using System.Linq;
using Microsoft.VisualStudio.TestTools.UnitTesting;
using RaytracerLib.MathLib;
using TomanuExtensions;
using RaytracerLib;
using System;
using RaytracerLib.MathLib.Polynomials;

namespace RaytracerTest
{
    [TestClass]
    public class SturmSequenceTest : TestBase
    {
        [TestMethod]
        public void SturmSequence_Test()
        {
            Polynomial p1 = new Polynomial(-1.0, 1.0) * 2;
            Polynomial p2 = new Polynomial(-2.0, 1.0) * -0.3;
            Polynomial p3 = new Polynomial(-3.0, 1.0) * -3;
            Polynomial p4 = new Polynomial(-4.0, 1.0) * 0.2;

            {
                Polynomial p = new Polynomial(4);
                SturmSequence sturm = new SturmSequence(p);

                Assert.IsTrue(sturm.Length == 1);
                Assert.IsTrue(sturm[0] == p);
                Assert.AreEqual(sturm.NumberOfSignChangesAtNegativeInfinity(), 0);
                Assert.AreEqual(sturm.NumberOfSignChangesAtPositiveInfinity(), 0);
            }

            {
                Polynomial p = p1;
                SturmSequence sturm = new SturmSequence(p);

                Assert.IsTrue(sturm.Length == 2);
                Assert.IsTrue(sturm[0] == p);
                Assert.IsTrue(sturm[1] == new Polynomial(2));
                Assert.AreEqual(sturm.NumberOfSignChangesAtNegativeInfinity(), 1);
                Assert.AreEqual(sturm.NumberOfSignChangesAtPositiveInfinity(), 0);
                Assert.AreEqual(sturm.NumberOfSignChangesAt(0.5), 1);
                Assert.AreEqual(sturm.NumberOfSignChangesAt(1.5), 0);
            }

            {

                Polynomial p = p1 * p1;
                SturmSequence sturm = new SturmSequence(p);

                Assert.IsTrue(sturm.Length == 3);
                Assert.IsTrue(sturm[0] == p);
                Assert.IsTrue(sturm[1] == new Polynomial(-8, 8));
                Assert.IsTrue(sturm[2] == new Polynomial(0));
                Assert.AreEqual(sturm.NumberOfSignChangesAtNegativeInfinity(), 1);
                Assert.AreEqual(sturm.NumberOfSignChangesAtPositiveInfinity(), 0);
                Assert.AreEqual(sturm.NumberOfSignChangesAt(0.5), 1);
                Assert.AreEqual(sturm.NumberOfSignChangesAt(1.5), 0);
            }

            {
                Polynomial p = p1 * p2;
                SturmSequence sturm = new SturmSequence(p);

                Assert.IsTrue(sturm.Length == 3);
                Assert.IsTrue(sturm[0] == p);
                Assert.IsTrue(sturm[1].IsAlmostRelativeEquals(new Polynomial(1.8, -1.2)));
                Assert.IsTrue(sturm[2].IsAlmostRelativeEquals(new Polynomial(-0.15)));
                Assert.AreEqual(sturm.NumberOfSignChangesAtNegativeInfinity(), 2);
                Assert.AreEqual(sturm.NumberOfSignChangesAtPositiveInfinity(), 0);
                Assert.AreEqual(sturm.NumberOfSignChangesAt(0.5), 2);
                Assert.AreEqual(sturm.NumberOfSignChangesAt(1.5), 1);
                Assert.AreEqual(sturm.NumberOfSignChangesAt(2.5), 0);
            }

            {
                Polynomial p = p1 * p2 * p3;
                SturmSequence sturm = new SturmSequence(p);

                Assert.IsTrue(sturm.Length == 4);
                Assert.IsTrue(sturm[0] == p);
                Assert.IsTrue(sturm[1].IsAlmostRelativeEquals(new Polynomial(19.8, -21.6, 5.4)));
                Assert.IsTrue(sturm[2].IsAlmostRelativeEquals(new Polynomial(-2.4, 1.2)));
                Assert.IsTrue(sturm[3].IsAlmostRelativeEquals(new Polynomial(1.8)));
                Assert.AreEqual(sturm.NumberOfSignChangesAtNegativeInfinity(), 3);
                Assert.AreEqual(sturm.NumberOfSignChangesAtPositiveInfinity(), 0);
                Assert.AreEqual(sturm.NumberOfSignChangesAt(0.5), 3);
                Assert.AreEqual(sturm.NumberOfSignChangesAt(1.5), 2);
                Assert.AreEqual(sturm.NumberOfSignChangesAt(2.5), 1);
                Assert.AreEqual(sturm.NumberOfSignChangesAt(3.5), 0);
            }

            {
                Polynomial p = p1 * p2 * p2;
                SturmSequence sturm = new SturmSequence(p);

                Assert.IsTrue(sturm.Length == 4);
                Assert.IsTrue(sturm[0] == p);
                Assert.IsTrue(sturm[1].IsAlmostRelativeEquals(new Polynomial(1.44, -1.8, 0.54)));
                Assert.IsTrue(sturm[2].IsAlmostRelativeEquals(new Polynomial(-0.08, 0.04)));
                Assert.IsTrue(sturm[3].IsAlmostRelativeEquals(new Polynomial(0)));
                Assert.AreEqual(sturm.NumberOfSignChangesAtNegativeInfinity(), 2);
                Assert.AreEqual(sturm.NumberOfSignChangesAtPositiveInfinity(), 1);
                Assert.AreEqual(sturm.NumberOfSignChangesAt(0.5), 2);
                Assert.AreEqual(sturm.NumberOfSignChangesAt(1.5), 1);
                Assert.AreEqual(sturm.NumberOfSignChangesAt(2.5), 1);
            }

            {
                Polynomial p = p2 * p2 * p2;
                SturmSequence sturm = new SturmSequence(p2 * p2 * p2);

                Assert.IsTrue(sturm.Length == 3);
                Assert.IsTrue(sturm[0] == p);
                Assert.IsTrue(sturm[1].IsAlmostRelativeEquals(new Polynomial(-0.324, 0.324, -0.081)));
                Assert.IsTrue(sturm[2].IsAlmostRelativeEquals(new Polynomial(0)));
                Assert.AreEqual(sturm.NumberOfSignChangesAtNegativeInfinity(), 2);
                Assert.AreEqual(sturm.NumberOfSignChangesAtPositiveInfinity(), 1);
                Assert.AreEqual(sturm.NumberOfSignChangesAt(1.5), 2);
                Assert.AreEqual(sturm.NumberOfSignChangesAt(2.5), 1);
            }

            {
                Polynomial p = p1 * p2 * p3 * p4;
                SturmSequence sturm = new SturmSequence(p);

                Assert.IsTrue(sturm.Length == 5);
                Assert.IsTrue(sturm[0] == p);
                Assert.IsTrue(sturm[1].IsAlmostRelativeEquals(new Polynomial(-18, 25.2, -10.8, 1.44)));
                Assert.IsTrue(sturm[2].IsAlmostRelativeEquals(new Polynomial(2.61, -2.25, 0.45)));
                Assert.IsTrue(sturm[3].IsAlmostRelativeEquals(new Polynomial(-2.88, 1.152)));
                Assert.IsTrue(sturm[4].IsAlmostRelativeEquals(new Polynomial(0.2025)));
                Assert.AreEqual(sturm.NumberOfSignChangesAtNegativeInfinity(), 4);
                Assert.AreEqual(sturm.NumberOfSignChangesAtPositiveInfinity(), 0);
                Assert.AreEqual(sturm.NumberOfSignChangesAt(0.5), 4);
                Assert.AreEqual(sturm.NumberOfSignChangesAt(1.5), 3);
                Assert.AreEqual(sturm.NumberOfSignChangesAt(2.5), 2);
                Assert.AreEqual(sturm.NumberOfSignChangesAt(3.5), 1);
                Assert.AreEqual(sturm.NumberOfSignChangesAt(4.5), 0);
            }

            {
                Polynomial p = p1 * p2 * p3 * p3;
                SturmSequence sturm = new SturmSequence(p);

                Assert.IsTrue(sturm.Length == 5);
                Assert.IsTrue(sturm[0] == p);
                Assert.IsTrue(sturm[1].IsAlmostRelativeEquals(new Polynomial(210.6, -313.2, 145.8, -21.6)));
                Assert.IsTrue(sturm[2].IsAlmostRelativeEquals(new Polynomial(-21.2625, 18.225, -3.7125)));
                Assert.IsTrue(sturm[3].IsAlmostRelativeEquals(new Polynomial(17.1371900826464, -5.71239669421584)));
                Assert.IsTrue(sturm[4].IsAlmostRelativeEquals(new Polynomial(0)));
                Assert.AreEqual(sturm.NumberOfSignChangesAtNegativeInfinity(), 4);
                Assert.AreEqual(sturm.NumberOfSignChangesAtPositiveInfinity(), 0);
                Assert.AreEqual(sturm.NumberOfSignChangesAt(0.5), 4);
                Assert.AreEqual(sturm.NumberOfSignChangesAt(1.5), 3);
                Assert.AreEqual(sturm.NumberOfSignChangesAt(2.5), 2);
                Assert.AreEqual(sturm.NumberOfSignChangesAt(3.5), 0);
            }

            {


                Polynomial p = p1 * p2 * p2 * p2;
                SturmSequence sturm = new SturmSequence(p);

                Assert.IsTrue(sturm.Length == 5);
                Assert.IsTrue(sturm[0] == p);
                Assert.IsTrue(sturm[1].IsAlmostRelativeEquals(new Polynomial(1.08, -1.944, 1.134, -0.216)));
                Assert.IsTrue(sturm[2].IsAlmostRelativeEquals(new Polynomial(-0.0405, 0.0405, -0.010125)));
                Assert.IsTrue(sturm[3].IsAlmostRelativeEquals(new Polynomial(-2.19824158875781E-14, 1.33226762955019E-14)));
                Assert.IsTrue(sturm[4].IsAlmostRelativeEquals(new Polynomial(0.00124031249999985)));
                Assert.AreEqual(sturm.NumberOfSignChangesAtNegativeInfinity(), 3);
                Assert.AreEqual(sturm.NumberOfSignChangesAtPositiveInfinity(), 1);
                Assert.AreEqual(sturm.NumberOfSignChangesAt(0.5), 3);
                Assert.AreEqual(sturm.NumberOfSignChangesAt(1.5), 2);
                Assert.AreEqual(sturm.NumberOfSignChangesAt(2.5), 1);
            }

            {
                Polynomial p = p1 * p1 * p1 * p1;
                SturmSequence sturm = new SturmSequence(p);

                Assert.IsTrue(sturm.Length == 3);
                Assert.IsTrue(sturm[0] == p);
                Assert.IsTrue(sturm[1].IsAlmostRelativeEquals(new Polynomial(-64, 192, -192, 64)));
                Assert.IsTrue(sturm[2].IsAlmostRelativeEquals(new Polynomial(0)));
                Assert.AreEqual(sturm.NumberOfSignChangesAtNegativeInfinity(), 1);
                Assert.AreEqual(sturm.NumberOfSignChangesAtPositiveInfinity(), 0);
                Assert.AreEqual(sturm.NumberOfSignChangesAt(0.5), 1);
                Assert.AreEqual(sturm.NumberOfSignChangesAt(1.5), 0);
            }

            {
                Polynomial p = p1 * p1 * p2 * p2;
                SturmSequence sturm = new SturmSequence(p);

                Assert.IsTrue(sturm.Length == 5);
                Assert.IsTrue(sturm[0] == p);
                Assert.IsTrue(sturm[1].IsAlmostRelativeEquals(new Polynomial(-4.32, 9.36, -6.48, 1.44)));
                Assert.IsTrue(sturm[2].IsAlmostRelativeEquals(new Polynomial(0.18, -0.27, 0.09)));
                Assert.IsTrue(sturm[3].IsAlmostRelativeEquals(new Polynomial(-1.33226762955019E-14, 1.4210854715202E-14)));
                Assert.IsTrue(sturm[4].IsAlmostRelativeEquals(new Polynomial(-0.00597656249999998)));
                Assert.AreEqual(sturm.NumberOfSignChangesAtNegativeInfinity(), 3);
                Assert.AreEqual(sturm.NumberOfSignChangesAtPositiveInfinity(), 1);
                Assert.AreEqual(sturm.NumberOfSignChangesAt(0.5), 3);
                Assert.AreEqual(sturm.NumberOfSignChangesAt(1.5), 3);
                Assert.AreEqual(sturm.NumberOfSignChangesAt(2.5), 1);
            }

            // Single roots.

            {
                Polynomial p = p1 * p1;
                Polynomial s = SturmSequence.SingleRoots(p).Normalize();
                Assert.IsTrue(s.IsAlmostRelativeEquals(p1.Normalize()));
            }

            {
                Polynomial p = p1 * p1 * p1;
                Polynomial s = SturmSequence.SingleRoots(p).Normalize();
                Assert.IsTrue(s.IsAlmostRelativeEquals(p1.Normalize()));
            }

            {
                Polynomial p = p1 * p1 * p1 * p1;
                Polynomial s = SturmSequence.SingleRoots(p).Normalize();
                Assert.IsTrue(s.IsAlmostRelativeEquals(p1.Normalize()));
            }

            {
                Polynomial p = p1 * p1 * p1 * p2;
                Polynomial s = SturmSequence.SingleRoots(p).Normalize();
                Assert.IsTrue(s.IsAlmostRelativeEquals((p1 * p2).Normalize()));
            }

            {
                Polynomial p = p1 * p1 * p2 * p2;
                Polynomial s = SturmSequence.SingleRoots(p).Normalize();
                Assert.IsTrue(s.IsAlmostRelativeEquals((p1 * p2).Normalize()));
            }

            {
                Polynomial p = p1 * p1 * p2 * p3;
                Polynomial s = SturmSequence.SingleRoots(p).Normalize();
                Assert.IsTrue(s.IsAlmostRelativeEquals((p1 * p2 * p3).Normalize()));
            }

            {
                Polynomial p = p1 * p1 * p2 * p3 * p4;
                Polynomial s = SturmSequence.SingleRoots(p).Normalize();
                Assert.IsTrue(s.IsAlmostRelativeEquals((p1 * p2 * p3 * p4).Normalize()));
            }

            {
                Polynomial p = p1 * p1 * p2 * p2 * p3 * p4;
                Polynomial s = SturmSequence.SingleRoots(p).Normalize();
                Assert.IsTrue(s.IsAlmostRelativeEquals((p1 * p2 * p3 * p4).Normalize()));
            }

            {
                Polynomial p = p1 * p1 * p2 * p2 * p3 * p3 * p4;
                Polynomial s = SturmSequence.SingleRoots(p, 1e-11).Normalize();
                Assert.IsTrue(s.IsAlmostRelativeEquals((p1 * p2 * p3 * p4).Normalize()));
            }

            {
                Polynomial p = p1 * p1 * p2 * p2 * p3 * p3 * p4 * p4;
                Polynomial s = SturmSequence.SingleRoots(p, 1e-9).Normalize();
                Assert.IsTrue(s.IsAlmostRelativeEquals((p1 * p2 * p3 * p4).Normalize(), 1e-10));
            }
        }
    }
}